Multiple Numerical Assembly 2

Sections

Abstract

Previously, I sent around some images comparing the alpha diversity results of each run. This met with general approval as well as questions including gamma diversity, the burstiness of arrivals (co-establishment), the burstiness of extinctions (co-extinctions), and comparisons with purely neutral models of arrival and extinction.

In this document, we continue to attempt to address Chris’s question of how increased spatial mobility changes biodiversity measures. We begin by loading up each data set referring to different neutral (arrival-extinction) and spatial (dispersal speed) conditions. We can then extract our diversity metrics for each system and each site/ecosystem/environment within the system and compare them.

Notes

• While Chris prefers abundance/richness as a dichotomy, it might be helpful for publication to consider more than that.

TODO

• Measure the impact of dispersal on invadability. This preferably identifies both the height of the new barrier as well as precisely determining the effect. This might include the asymmetry at the ends of the line spatial arrangement (as compared to a ring spatial arrangement). ** I anticipate that this is precisely the effect of the extra term. To show it then, analytically we can check the per-capita rate of growth with the additional term, which will act as a penalty. Empirically, we look at the invasions that we expect to have succeeded but then died out when space was introduced and we should be able to see a good match between its predicted per-capita growth rate and when/if it died out within a given system. This breaks down if it can survive in an adjacent system, for what it is worth.
• Evenness by trophic level. This might help us deduce how neutral the trophic levels are, but requires analysing how the trophic network changes over time (possibly in contrast to the theoretical food web). Furthermore, this is most advantageous to investigate as we change the parameters of the simulation to include more trophic levels. This would also allow us to color code abundance and examine richness of trophic levels (both number of levels and richness within levels). ** We should make a point of running another layer of consumers (and producers?) as well so as to be able to better observe the resultant diversity. I do note that, contrary to previous expecations, we are running at about \(30\%\) of the pool. As I am only using 1 pool at this point, there might be a high amount of variation. ** The parent point here might require moving onto Viking. The deSolve package does have functions that return lists, where the first element is a vector of derivatives of \(y\) with respect to \(t\). The ‘’next’’ elements seem like they could be things like the trophic structure and trophic levels.
• We should consider the “run to stability”, then variable change. This requires another set of (2 sets of) systems where we run the system as before to stability (everything looks stable at least!) and then make a jump in one (or two? space and arrival or space and extinction rate?) of the parameters.
• From previous chats, there is also the prospect of comparing ‘’pure’’ neutral dynamics. In terms of simulations, a pool of the same size, but all producers would seem to do it. In terms of analytics, there should be a good neutral theory solution.
• Another point of interest is how bursty is the distribution compared to the exponential null model? This in effect studies the impact of the dynamics. If they are far burstier than we expect/neutrality would dictate, than we can deduce there are multiple related secondary extinctions. Otherwise, the system would appear to be merely neutral dynamics and primary extinctions.

Future

• Once the strictly deterministic dynamics are sufficiently cleared out, Chris (and Susan) would like to swap to Gillespie style / jump process / stochiometry dynamics.

Functions

library(dplyr)        # Data Manipulation
library(tidyr)        # Data Pivotting
library(ggplot2)      # 2-D Plot
library(plotly)       # 3-D Plot
library(ggfortify)    # used for biplots of PCAs
library(vegan)        # Ecological analysis mega-package
library(htmltools)    # Programmatic chunks
library(fitdistrplus) # Advanced fitting over MASS.
library(poweRlaw)     # Fitting heavy-tailed distributions.
library(RMTRCode2)    # Personal package.
# https://stackoverflow.com/a/7172832
ifrm <- function(obj, env = globalenv()) {
  obj <- deparse(substitute(obj))
  if(exists(obj, envir = env)) {
    rm(list = obj, envir = env)
  }
}
log0 <- function(x, base = exp(1)) {
  xneq0 <- x != 0 & !is.na(x)
  x[xneq0] <- log(x[xneq0], base = base)
  return(x)
}

Data

As of the writing of this script, all of the data is stored in the same location as this script. The parameters used to generate the interaction matrices, pools, and events is not stored with the attempts, which is something I should probably fix in the future. (I could store them in the ellipsis argument or as a separate object.)

by_for_thinning <- 10 # time steps
divide_time_by <- 1E4 # time units
burn_in <- 1E4 # time units
load_safe_thin <- function(fname, bythin, divtime, burn) {
  loaded <- tryCatch({load(fname)}, 
                     error = function(e) {
                       print(fname)
                       print(e)
                       return(NA)
                     })
  if (is.na(loaded)) {
    return(NA)
  } else {
    loaded <- get(loaded)
    loaded$Abundance <- loaded$Abundance[seq(from = 1,
                                             to = nrow(loaded$Abundance),
                                             by = bythin), ]
    
    toEliminate <- loaded$Abundance[, -1] < loaded$Parameters$EliminationThreshold & loaded$Abundance[, -1] > 0
    loaded$Abundance[, -1][toEliminate] <- 0
    
    loaded$Abundance <- loaded$Abundance[
      loaded$Abundance[, 1] > burn,
    ]
    
    loaded$Abundance[, 1] <- loaded$Abundance[, 1] / divtime
    return(loaded)
  }
}
# All .RData
files_dat <- dir(pattern = ".RData$")
# Remove PoolMats
files_dat <- files_dat[
  !grepl(x = files_dat, 
         pattern = "PoolMats", 
         fixed = TRUE)
  ]
# Technically overkill, but prevents unintentional loads.
# Break into two separate runs to load only intended.
# Process "MNA-FirstAttempt#####-Result-Env10-####.RData"
files_dat_FA <- files_dat[
  grepl(x = files_dat, 
        pattern = "FirstAttempt", 
        fixed = TRUE)
  ]
Results <- sapply(
  files_dat_FA,
  load_safe_thin, 
  bythin = by_for_thinning,
  divtime = divide_time_by, 
  burn = burn_in,
  simplify = FALSE, USE.NAMES = TRUE
)
# Process "MNA-Dist#####-Ext###-Env10-####.RData"
files_dat_Ext <- files_dat[
  grepl(x = files_dat, 
        pattern = "Dist", 
        fixed = TRUE)
  ]
Results <- c(Results, sapply(
  files_dat_Ext,
  load_safe_thin, 
  bythin = by_for_thinning, 
  divtime = divide_time_by, 
  burn = burn_in,
  simplify = FALSE, USE.NAMES = TRUE
))

Pools and Matrices

load_safe <- function(fname) {
  loaded <- tryCatch({load(fname)}, 
                     error = function(e) {
                       print(fname)
                       print(e)
                       return(NA)
                     })
  if (all(is.na(loaded))) {
    return(NA)
  } else {
    return(sapply(loaded, get, 
                  envir = sys.frame(sys.parent(0)), 
                  simplify = FALSE, USE.NAMES = TRUE))
  }
}
# All .RData
files_dat_PM <- dir(pattern = ".RData$")
# Remove PoolMats
files_dat_PM <- files_dat_PM[
  grepl(x = files_dat_PM, 
        pattern = "PoolMats", 
        fixed = TRUE)
  ]
PoolsMats <- sapply(
  files_dat_PM,
  load_safe, 
  simplify = FALSE, USE.NAMES = TRUE
)

Trophic Functions

EliminiationThreshold <- unique(sapply(
  Results, 
  function(lst) {lst$Parameters$EliminationThreshold}
))
stopifnot(length(EliminiationThreshold) == 1)
NumEnvironments <- unique(sapply(
  Results, 
  function(lst) {lst$NumEnvironments}
))
stopifnot(length(NumEnvironments) == 1)
TrophicFunctions <- sapply(
  PoolsMats,
  function(PM, NE, ET) {
    RMTRCode2::CalculateTrophicStructure(
      Pool = PM$Pool,
      NumEnvironments = NE,
      InteractionMatrices = PM$InteractionMatrices,
      EliminationThreshold = ET
    )
  },
  NE = NumEnvironments,
  ET = EliminiationThreshold
)

#TODO Fix warning (binding character and factor vector, coercing into character vector in bind_rows)
TrophicAnalyses <- list()
Result = Results
Nm = names(Results)
TrophFN = TrophicFunctions

for (i in seq_along(Results)) {
  print(i); print(Nm[i])
  # Identify Appropriate Function.
  Unusual1 <- grepl(pattern = "FirstAttempt",
                    x = Nm[i], fixed = TRUE)
  NormalKey <- strsplit(Nm[i], split = '-')[[1]][3]
  linkedFN <- if(Unusual1) {
    TrophFN[grepl(x = names(TrophFN),
                  pattern = "FirstAttempt", fixed = TRUE)][[1]]
  } else {
    TrophFN[grepl(x = names(TrophFN),
                  pattern = NormalKey, fixed = TRUE)][[1]]
  }
  
  # Apply to each row.
  
  TrophicAnalyses[[i]] <- apply(Result[[i]]$Abundance[, -1],
                                MARGIN = 1,
                                FUN = linkedFN)
}

# TrophicAnalyses <- lapply(
#   seq_along(Results),
#   function(i, Result, Nm, TrophFN) {
#     # Identify Appropriate Function.
#     Unusual1 <- grepl(pattern = "FirstAttempt",
#                       x = Nm[i], fixed = TRUE)
#     NormalKey <- strsplit(Nm[i], split = '-')[[1]][3]
#     linkedFN <- if(Unusual1) {
#       TrophFN[grepl(x = names(TrophFN),
#                      pattern = "FirstAttempt", fixed = TRUE)][[1]]
#     } else {
#       TrophFN[grepl(x = names(TrophFN), 
#                      pattern = NormalKey, fixed = TRUE)][[1]]
#     }
#     
#     # Apply to each row.
#     
#     return(apply(Result[[i]]$Abundance[, -1],
#                  MARGIN = 1,
#                  FUN = linkedFN)
#     )
#   },
#   Result = Results,
#   Nm = names(Results),
#   TrophFN = TrophicFunctions
# )

names(TrophicAnalyses) <- names(Results)

Diversity

Preparation

# Borrowing code from FirstAttempt-Doc-Analysis.Rmd
Calculate_Diversity <- function(result) {
  Diversity <- lapply(
    1:result$NumEnvironments,
    function(i, abund, numSpecies) {
      time <- abund[, 1]
      env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
      richness <- rowSums(env != 0)
      abundSum <- rowSums(env)
      #NOTE: THIS CAN YIELD NAN'S (0/0).
      # THIS IS NOT NECESSARILY A PROBLEM.
      # IT MIGHT BE WORTH IT JUST TO USE 0 OR
      # TO CATCH IT EXPLICITLY AND REPLACE WITH NAN.
      entropy <- env / abundSum
      entropy <- - apply(
        entropy, MARGIN = 1,
        FUN = function(x) {
          sum(x * log0(x))
        })
      species <- apply(
        env, MARGIN = 1,
        FUN = function(x) {
          toString(which(x > 0))
        }
      )
      evenness <- entropy / log(richness)
      data.frame(Time = time, 
                 Richness = richness, 
                 Entropy = entropy,
                 Evenness = evenness,
                 Species = species,
                 Environment = i,
                 stringsAsFactors = FALSE)
    },
    abund = result$Abundance,
    numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
  )
  
  
  Diversity <- dplyr::bind_rows(Diversity)
  Diversity_alpha <- Diversity
  # Diversity_alpha <- Diversity_alpha %>% dplyr::mutate(
  #   Evenness = Entropy / log(Richness)
  # )
  
  # Modify to do the gamma bits right here.
  Diversity_gamma <- Diversity %>% dplyr::group_by(
    Time
  ) %>% dplyr::summarise(
    Mean = mean(Richness),
    SpeciesTotal = toString(sort(unique(unlist(strsplit(paste(
      Species, collapse = ", "), split = ", ", fixed = TRUE))))),
    Gamma = unlist(lapply(strsplit(
      SpeciesTotal, split = ", ", fixed = TRUE), function(x) length(x[x!=""]) ))
  ) %>% tidyr::pivot_longer(
    cols = c(Mean, Gamma), 
    names_to = "Aggregation",
    values_to = "Richness"
  )
  
  # Combine the two types of results
  Diversity_alpha <- Diversity_alpha %>% dplyr::select(
    -Species
  ) %>% tidyr::pivot_longer(
    cols = c(Richness, Entropy, Evenness),
    names_to = "Measurement",
    values_to = "Value"
  ) %>% dplyr::mutate(
    Environment = as.character(Environment)
  )
  
  Diversity_gamma <- Diversity_gamma %>% dplyr::select(
    -SpeciesTotal
  ) %>% dplyr::rename(
    Environment = Aggregation,
    Value = Richness
  ) %>% dplyr::mutate(
    Measurement = "Richness"
  )
  
  Diversity_beta <- Diversity_alpha %>% dplyr::filter(
    Measurement == "Richness"
  ) %>% dplyr::select(
    -Measurement
  ) %>% dplyr::left_join(
    y = Diversity_gamma %>% dplyr::filter(
      Measurement == "Richness", Environment == "Gamma"
    ) %>% dplyr::select(
      -Measurement, -Environment
    ),
    by = "Time",
    suffix = c("_Alpha", "_Gamma")
    # ) %>% dplyr::group_by(
    #   Time
  ) %>% dplyr::mutate(
    BetaSpeciesMissing = Value_Gamma - Value_Alpha,
    BetaSpeciesPercentage = Value_Alpha/Value_Gamma
  ) %>% dplyr::select(
    -Value_Gamma, -Value_Alpha
  ) %>% tidyr::pivot_longer(
    names_to = "Measurement",
    values_to = "Value",
    cols = c(BetaSpeciesMissing, BetaSpeciesPercentage)
    # ) %>% dplyr::ungroup(
  )
  
  #print(c(colnames(Diversity_alpha), colnames(Diversity_beta), colnames(Diversity_gamma)))
  Diversity <- rbind(
    Diversity_alpha,
    Diversity_beta,
    Diversity_gamma
  )
  
  return(Diversity)
}
Diversity_jaccard_space <- function(result) {
  apply(
    result$Abundance,
    MARGIN = 1, # Rows
    function(row, envs) {
      time <- row[1]
      dists <- vegan::vegdist(
        method = "jaccard", 
        x = matrix(row[-1] > 0, nrow = envs, byrow = TRUE)
      )
      
      dataf <- expand.grid(
        Env1 = 1:envs,
        Env2 = 1:envs
      ) %>% dplyr::filter(
        Env1 < Env2
      ) %>% dplyr::mutate(
        Time = time,
        Jaccard = dists
      )
      
      return(dataf)
    },
    envs = result$NumEnvironments
  )
}
Diversity_jaccard_time <- function(result, subsample = 100) {
  # Break into environments, then apply it to the time series.
  patches <- lapply(
    1:result$NumEnvironments, function(i, abund, envs, spec) {
      abund <- abund[seq(from = 1, to = nrow(abund), by = subsample), ]
      times <- abund[, 1]
      patch <- abund[, 1 + 1:spec + spec * (i - 1)]
      
      dists <- vegan::vegdist(
        method = "jaccard", 
        x = patch > 0
      )
      
      dataf <- expand.grid(
        Time1 = times,
        Time2 = times
      ) %>% dplyr::filter(
        Time1 < Time2
      ) %>% dplyr::mutate(
        Environment = i,
        Jaccard = dists
      )
      
      return(dataf)
    }, abund = result$Abundance, envs = result$NumEnvironments,
    spec = (ncol(result$Abundance) - 1) / result$NumEnvironments
  )
}

Calculate_Species <- function(result, bintimes = FALSE) {
  SpeciesPerEnvironment <- lapply(
    1:result$NumEnvironments,
    function(i, abund, numSpecies) {
      time <- abund[, 1]
      env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
      # Need to retrieve Position and Value
      species <- apply(
        cbind(time, env), MARGIN = 1,
        FUN = function(x) {
          time <- x[1]
          dat <- x[-1]
          if (any(dat > 0)) {
            positions <- (which(dat > 0))
            values <- dat[positions]
            data.frame(
              Time = time,
              Species = positions,
              Abundance = values,
              row.names = NULL
            )
          } else {NULL}
          # Returns as list
        }
      )
      return(
        dplyr::bind_rows(species) %>% dplyr::mutate(
          Environment = i
        )
      )
    },
    abund = result$Abundance,
    numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
  )
  
  if (bintimes) {
    # Should equalise time steps.
    SpeciesPerEnvironment <- lapply(
      SpeciesPerEnvironment, function(SPE) {
        SPE %>% dplyr::mutate(
          TimeFloor = floor(Time*10)/10
        ) %>% dplyr::group_by(
          TimeFloor, Species, Environment
        ) %>% dplyr::summarise(
          Abundance = median(Abundance, na.rm = TRUE)
        )
      })
  }
  
  return(dplyr::bind_rows(SpeciesPerEnvironment))
}

# Note that if a file fails to load, we might have NA instead of a result to work with.
Diversity <- sapply(
  USE.NAMES = TRUE, simplify = FALSE,
  Results, function(result) {
    if (length(result) == 1 && is.na(result)) {
      # Problem case.
      return(NA)
    }
    # print(paste("Calculating", Sys.time()))
    
    # Calculate the diversity.
    # We will need to extract the system properties from
    # the file names which we carry through using sapply.
    return(Calculate_Diversity(result))
  }
)

# Expect warnings since we have all 0 rows on occasion.
JaccardSpace <- sapply(
  USE.NAMES = TRUE, simplify = FALSE,
  Results, function(result) {
    if (length(result) == 1 && is.na(result)) {
      # Problem case.
      return(NA)
    }
    # print(paste("Calculating", Sys.time()))
    
    # Calculate the diversity.
    # We will need to extract the system properties from
    # the file names which we carry through using sapply.
    return(Diversity_jaccard_space(result))
  }
)

# Expect warnings since we have all 0 rows on occasion.
JaccardTime <- sapply(
  USE.NAMES = TRUE, simplify = FALSE,
  Results, function(result) {
    if (length(result) == 1 && is.na(result)) {
      # Problem case.
      return(NA)
    }
    print(paste("Calculating", Sys.time()))
    
    # Calculate the diversity.
    # We will need to extract the system properties from
    # the file names which we carry through using sapply.
    return(Diversity_jaccard_time(result))
    # Too many time points for calculations.
    # Needs to have log(length(times)^2, base = 2) < 31.
    # => length(times) < 46340 or so.
  }
)

[1] "Calculating 2021-11-14 18:47:55"
[1] "Calculating 2021-11-14 18:47:55"
[1] "Calculating 2021-11-14 18:47:55"
[1] "Calculating 2021-11-14 18:47:55"
[1] "Calculating 2021-11-14 18:47:56"
[1] "Calculating 2021-11-14 18:47:56"
[1] "Calculating 2021-11-14 18:47:56"
[1] "Calculating 2021-11-14 18:47:56"
[1] "Calculating 2021-11-14 18:47:56"
[1] "Calculating 2021-11-14 18:47:57"
[1] "Calculating 2021-11-14 18:47:57"
[1] "Calculating 2021-11-14 18:47:57"
[1] "Calculating 2021-11-14 18:47:57"
[1] "Calculating 2021-11-14 18:47:58"
[1] "Calculating 2021-11-14 18:47:59"
[1] "Calculating 2021-11-14 18:47:59"
[1] "Calculating 2021-11-14 18:47:59"
[1] "Calculating 2021-11-14 18:47:59"
[1] "Calculating 2021-11-14 18:48:00"
[1] "Calculating 2021-11-14 18:48:00"
[1] "Calculating 2021-11-14 18:48:01"
[1] "Calculating 2021-11-14 18:48:01"
[1] "Calculating 2021-11-14 18:48:01"
[1] "Calculating 2021-11-14 18:48:01"
[1] "Calculating 2021-11-14 18:48:02"
[1] "Calculating 2021-11-14 18:48:02"

SpeciesPresence <-  sapply(
  USE.NAMES = TRUE, simplify = FALSE,
  Results, function(result) {
    if (length(result) == 1 && is.na(result)) {
      # Problem case.
      return(NA)
    }
    # print(paste("Calculating", Sys.time()))
    
    # Calculate the diversity.
    # We will need to extract the system properties from
    # the file names which we carry through using sapply.
    return(Calculate_Species(result))
  }
)

Properties <- strsplit(names(Diversity), '-', 
                       fixed = TRUE)
# 1st Chunk: Name, Discard
# 2nd Chunk: Iteration + Distance
# 3rd Chunk: Result or Extinction Rate (or Arrival?)
# 4th Chunk: Number of Environments
# 5th Chunk: Space Type + .RData
# Note the mix of Keyword and Location Structure (oops).
# Note also that this strsplit character is a bad decision and should be changed for next time. (D'oh.)
# (E.g. Dates DD-MM-YYYY, Decimals 1.35e-05.)
Properties <- data.frame(
  do.call(rbind, Properties),
  stringsAsFactors = FALSE
)
names(Properties)[1:5] <- c(
  "Name", "IterANDDist", "Modifier", "EnvNum", "SpaceAND.RData"
)
Properties$FullName <- names(Diversity)
# Capture the position between the text (first group)
# and the set of numbers (somehow without the +).
# The \\K resets so that we do not capture any text.
patternString <- "((?>[a-zA-Z]+)(?=[0-9eE]))\\K"
# Split strings. Some of the trick will be to introduce
# a character to make the separation around. We use "_".
Properties <- Properties %>% dplyr::mutate(
  IterANDDist = gsub(pattern = patternString, 
                     replacement = "_", 
                     x = IterANDDist, perl = TRUE),
  Modifier = gsub(pattern = patternString, 
                  replacement = "_", 
                  x = Modifier, perl = TRUE),
  EnvNum = gsub(pattern = patternString, 
                replacement = "_", 
                x = EnvNum, perl = TRUE)
) %>% tidyr::separate(
  IterANDDist, into = c("Iter", "Distance"),
  sep = "[_]", fill = "right"
) %>% tidyr::separate(
  Modifier, into = c("Modifier", "ModIntensity"),
  sep = "[_]", fill = "right"
) %>% tidyr::separate(
  EnvNum, into = c("Env", "Environments"),
  sep = "[_]"
) %>% tidyr::separate(
  SpaceAND.RData, into = c("Space", ".RData"),
  sep = "[.]"
) %>% dplyr::select(
  -Name, -.RData, -Env
) %>% dplyr::mutate(
  Distance = dplyr::case_when(
    is.na(Distance) ~ "1e+00",
    TRUE ~ Distance
  )
)

Diversity <- lapply(1:length(Diversity),
                    function(i, df, nm) {
                      df[[i]] %>% mutate(
                        Simulation = nm[i]
                      )
                    }, 
                    df = Diversity,
                    nm = names(Diversity))

JaccardSpace <- lapply(1:length(JaccardSpace),
                    function(i, df, nm) {
                      df[[i]] %>% dplyr::bind_rows(
                        # Need to account for the by times...
                        # Generates attribute warnings.
                      ) %>% mutate(
                        Simulation = nm[i]
                      )
                    }, 
                    df = JaccardSpace,
                    nm = names(JaccardSpace))

JaccardTime <- lapply(1:length(JaccardTime),
                    function(i, df, nm) {
                      df[[i]] %>% dplyr::bind_rows(
                        # Need to account for the by envs...
                        # Generates attribute warnings.
                      ) %>% mutate(
                        Simulation = nm[i]
                      )
                    }, 
                    df = JaccardTime,
                    nm = names(JaccardTime))

SpeciesPresence <- lapply(1:length(SpeciesPresence),
                    function(i, df, nm) {
                      df[[i]] %>% mutate(
                        Simulation = nm[i]
                      )
                    }, 
                    df = SpeciesPresence,
                    nm = names(SpeciesPresence))

Diversity <- dplyr::left_join(
  dplyr::bind_rows(Diversity),
  Properties,
  by = c("Simulation" = "FullName")
)

JaccardSpace <- dplyr::left_join(
  dplyr::bind_rows(JaccardSpace),
  Properties,
  by = c("Simulation" = "FullName")
)

JaccardTime <- dplyr::left_join(
  dplyr::bind_rows(JaccardTime),
  Properties,
  by = c("Simulation" = "FullName")
)

SpeciesPresence <- dplyr::left_join(
  dplyr::bind_rows(SpeciesPresence),
  Properties,
  by = c("Simulation" = "FullName")
)

Diversity <- Diversity %>% dplyr::mutate(
  Distance = dplyr::case_when(
    is.na(Distance) ~ "1e+00",
    TRUE ~ Distance
  )
)

Alpha Richness

Overall

ggplot2::ggplot(
  Diversity %>% dplyr::filter(
    Measurement == "Richness",
    Environment != "Gamma",
    Space != "Ring",
    Environment != "Mean"
  ), 
  ggplot2::aes(
    x = Time,
    y = Value,
    color = factor(Environment),
    alpha = ifelse(Environment %in% c("3", "7", "Mean"), 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::filter(
    Measurement == "Richness",
    Environment != "Gamma",
    Space != "Ring",
    Environment == "Mean"
  ), 
  color = "black"
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
) + ggplot2::labs(
  y = "Richness",
  x = paste0("Time, ", divide_time_by, " units"),
  title = "Overview of Alpha Richness over Time by System Properties",
  caption = "Each row has the same arrival and extinction events."
) + ggplot2::facet_grid(
  Modifier + ModIntensity ~ Space + Distance
)

# ggplot2::ggsave(overallrich + ggplot2::coord_cartesian(ylim = c(0, 25)), filename = "MNA-AlphaRichness-Overview.pdf", dpi = "retina", width = 11, height = 8)

Plot_Richness <- function(df) {
  tempname <- paste(unique(df$Simulation), collapse = " ")
  
  temp <- ggplot2::ggplot(
    df %>% dplyr::filter(
      Measurement == "Richness",
      Environment != "Gamma",
      Environment != "Mean"
    ), 
    ggplot2::aes(
      x = Time,
      y = Value,
      color = factor(Environment),
      alpha = ifelse(Environment %in% c("3", "7", "Mean"), 1, 0.3)
    )
  ) + ggplot2::geom_line(
  ) + ggplot2::geom_line(
    data = df %>% dplyr::filter(
      Measurement == "Richness",
      Environment != "Gamma",
      Environment == "Mean"
    ), 
    color = "black"
  ) + ggplot2::guides(
    alpha = "none"
  ) + ggplot2::scale_color_discrete(
    "Environment"
  ) + ggplot2::labs(
    title = tempname,
    x = paste0("Time, ", divide_time_by, " units"),
    y = "Richness"
  )
  
  return(temp)
}
Plots_Richness_alpha <- Diversity %>% dplyr::group_split(
  Simulation
) %>% purrr::map(
  Plot_Richness
)
# In the next chunk, we use
# www.r-bloggers.com/2020/07/programmatically-create-new-geadings-and-outputs-in-rmarkdown/

Sim: 1

Sim: 2

Sim: 3

Sim: 4

Sim: 5

Sim: 6

Sim: 7

Sim: 8

Sim: 9

Sim: 10

Sim: 11

Sim: 12

Sim: 13

Sim: 14

Sim: 15

Sim: 16

Sim: 17

Sim: 18

Sim: 19

Sim: 20

Sim: 21

Sim: 22

Sim: 23

Sim: 24

Sim: 25

Sim: 26

Evenness

Overall

ggplot2::ggplot(
  Diversity %>% dplyr::filter(
    Measurement == "Evenness",
    Environment != "Gamma",
    Space != "Ring",
    Environment != "Mean"
  ), 
  ggplot2::aes(
    x = Time,
    y = Value,
    color = factor(Environment),
    alpha = ifelse(Environment %in% c("3", "7", "Mean"), 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::filter(
    Measurement == "Evenness",
    Environment != "Gamma",
    Space != "Ring",
    Environment == "Mean"
  ), 
  color = "black"
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
) + ggplot2::labs(
  y = "Evenness",
  x = paste0("Time, ", divide_time_by, " units"),
  title = "Overview of Evenness over Time by System Properties",
  caption = "Each row has the same arrival and extinction events."
) + ggplot2::facet_grid(
  Modifier + ModIntensity ~ Space + Distance
)

# ggplot2::ggsave(overallrich + ggplot2::coord_cartesian(ylim = c(0, 25)), filename = "MNA-AlphaRichness-Overview.pdf", dpi = "retina", width = 11, height = 8)

Plot_Evenness <- function(df) {
  tempname <- paste(unique(df$Simulation), collapse = " ")
  
  temp <- ggplot2::ggplot(
    df %>% dplyr::filter(
      Measurement == "Evenness",
      Environment != "Gamma",
      Environment != "Mean"
    ), 
    ggplot2::aes(
      x = Time,
      y = Value,
      color = factor(Environment),
      alpha = ifelse(Environment %in% c("3", "7", "Mean"), 1, 0.3)
    )
  ) + ggplot2::geom_line(
  ) + ggplot2::geom_line(
    data = df %>% dplyr::filter(
      Measurement == "Evenness",
      Environment != "Gamma",
      Environment == "Mean"
    ), 
    color = "black"
  ) + ggplot2::guides(
    alpha = "none"
  ) + ggplot2::scale_color_discrete(
    "Environment"
  ) + ggplot2::labs(
    title = tempname,
    x = paste0("Time, ", divide_time_by, " units"),
    y = "Evenness"
  )
  
  return(temp)
}
Plots_Evenness <- Diversity %>% dplyr::group_split(
  Simulation
) %>% purrr::map(
  Plot_Evenness
)
# In the next chunk, we use
# www.r-bloggers.com/2020/07/programmatically-create-new-headings-and-outputs-in-rmarkdown/

Sim: 1

Sim: 2

Sim: 3

Sim: 4

Sim: 5

Sim: 6

Sim: 7

Sim: 8

Sim: 9

Sim: 10

Sim: 11

Sim: 12

Sim: 13

Sim: 14

Sim: 15

Sim: 16

Sim: 17

Sim: 18

Sim: 19

Sim: 20

Sim: 21

Sim: 22

Sim: 23

Sim: 24

Sim: 25

Sim: 26

Gamma Richness

Richness

ggplot2::ggplot(
  Diversity %>% dplyr::filter(
    Measurement == "Richness",
    Environment == "Gamma",
    Space != "Ring"
  ), 
  ggplot2::aes(
    x = Time,
    y = Value
  )
) + ggplot2::geom_line(
) + ggplot2::labs(
  y = "Richness",
  x = paste0("Time, ", divide_time_by, " units"),
  title = "Gamma Richness over Time by System Properties",
  caption = "Each row has the same arrival and extinction events."
) + ggplot2::facet_grid(
  Modifier + ModIntensity ~ Space + Distance
) + ggplot2::coord_cartesian(ylim = c(0, 45))

# ggplot2::ggsave(overallgamma, filename = "MNA-GammaRichness-Overview.pdf", dpi = "retina", width = 11, height = 8)

Beta Richness

\(\gamma - \alpha\)

ggplot2::ggplot(
  Diversity %>% dplyr::filter(
    Measurement == "BetaSpeciesMissing",
    Environment != "Gamma",
    Space != "Ring",
    Environment != "Mean"
  ), 
  ggplot2::aes(
    x = Time,
    y = Value,
    color = factor(Environment),
    alpha = ifelse(Environment %in% c("3", "7", "Mean"), 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
    data = Diversity %>% dplyr::filter(
        Measurement == "BetaSpeciesMissing",
        Environment != "Gamma",
        Space != "Ring",
        Environment != "Mean"
    ) %>% dplyr::group_by(
        Modifier, ModIntensity, Space, Distance, Time
    ) %>% dplyr::summarise(
        Environment = "Mean",
        Value = mean(Value, na.rm = TRUE)
    ), color = "black"
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
) + ggplot2::labs(
  y = "Absolute Species Turnover", # en.wikipedia.org/wiki/Beta_diversity
  x = paste0("Time, ", divide_time_by, " units"),
  title = "Overview of Gamma - Alpha Richness over Time by System Properties",
  caption = "Each row has the same arrival and extinction events."
) + ggplot2::facet_grid(
  Modifier + ModIntensity ~ Space + Distance
)

# ggplot2::ggsave(overallbetamiss, filename = "MNA-BetaMissing-Overview.pdf", dpi = "retina", width = 11, height = 8)

\(\alpha / \gamma\)

ggplot2::ggplot(
    Diversity %>% dplyr::filter(
        Measurement == "BetaSpeciesPercentage",
        Environment != "Gamma",
        Space != "Ring",
        Environment != "Mean"
    ), 
    ggplot2::aes(
        x = Time,
        y = Value,
        color = factor(Environment),
        alpha = ifelse(Environment %in% c("3", "7", "Mean"), 1, 0.3)
    )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
    data = Diversity %>% dplyr::filter(
        Measurement == "BetaSpeciesPercentage",
        Environment != "Gamma",
        Space != "Ring",
        Environment != "Mean"
    ) %>% dplyr::group_by(
        Modifier, ModIntensity, Space, Distance, Time
    ) %>% dplyr::summarise(
        Environment = "Mean",
        Value = mean(Value, na.rm = TRUE)
    ), color = "black"
) + ggplot2::guides(
    alpha = "none"
) + ggplot2::scale_color_discrete(
    "Environment"
) + ggplot2::labs(
    y = "Percentage Species Present", # en.wikipedia.org/wiki/Beta_diversity
    x = paste0("Time, ", divide_time_by, " units"),
    title = "Overview of Alpha/Gamma Richness over Time by System Properties",
    caption = "Each row has the same arrival and extinction events."
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
)

# ggplot2::ggsave(overallbetapercent, filename = "MNA-BetaPercent-Overview.pdf", dpi = "retina", width = 11, height = 8)

Alpha vs Gamma

Diversity_AG <- dplyr::left_join(
  Diversity %>% dplyr::filter(
    Measurement == "Richness",
    Space != "Ring",
    Environment != "Gamma"
  ), 
  Diversity %>% dplyr::filter(
    Measurement == "Richness",
    Space != "Ring",
    Environment == "Gamma"
  ) %>% dplyr::rename(
    Gamma = Value
  ) %>% dplyr::select(
    -Environment
  )
) %>% dplyr::mutate(
  Distance = ifelse(is.na(Distance), 1, Distance)
)

Joining, by = c("Time", "Measurement", "Simulation", "Iter", "Distance", "Modifier", "ModIntensity", "Environments", "Space")

Diversity_AG_Binned <- Diversity_AG %>% dplyr::mutate(
    TimeFloor = floor(Time * 10) / 10
) %>% dplyr::distinct(
    Modifier, ModIntensity, Space, Distance, # Simulation/Facets
    Environment, TimeFloor, # Grouping Bins
    Value, Gamma # Values. If either moves, the trajectory entered a new square.
)
# ggplot2::ggplot(
#   temp %>% dplyr::filter(
#     Modifier == "Result",
#     Space == "None",
#     Environment != "Mean"
#   ), ggplot2::aes(
#     x = Gamma,
#     y = Value
#   )
# ) + ggplot2::geom_bin2d(
# ) + ggplot2::scale_fill_viridis_c(
# ) + ggplot2::labs(
#   Title = "No Spatial Structure, Equal Extinction & Arrival Rates",
#   x = "Gamma Richness",
#   y = "Alpha Richness"
# )

ggplot2::ggplot(
    Diversity_AG %>% dplyr::filter(
        Environment != "Mean",
        Space != "Ring"
    ), ggplot2::aes(
        x = Value,
        y = Gamma
    )
) + ggplot2::geom_bin2d(
) + ggplot2::scale_fill_viridis_c(
    trans = "log10"
) + ggplot2::geom_point(
    data = Diversity_AG %>% dplyr::filter(
        Environment != "Mean",
        Space != "Ring"
    ) %>% dplyr::group_by(
        Modifier, ModIntensity, Space, Distance
    ) %>% dplyr::summarise(
        Gamma = mean(Gamma, na.rm = TRUE),
        Value = mean(Value, na.rm = TRUE)
    ), color = "red", size = 2, shape = 4
) + ggplot2::labs(
    x = "Alpha Richness",
    y = "Gamma Richness"
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::geom_abline(slope = 1, intercept = 0)

# ggplot2::ggsave(overallalphagamma, filename = "MNA-AlphaGamma-Overview.pdf", dpi = "retina", width = 11, height = 8)

ggplot2::ggplot(
    Diversity_AG_Binned %>% dplyr::filter(
        Environment != "Mean",
        Space != "Ring"
    ), ggplot2::aes(
        x = Value,
        y = Gamma
    )
) + ggplot2::geom_bin2d(
) + ggplot2::scale_fill_viridis_c(
    trans = "log10"
) + ggplot2::geom_point(
    data = Diversity_AG_Binned %>% dplyr::filter(
        Environment != "Mean",
        Space != "Ring"
    ) %>% dplyr::group_by(
        Modifier, ModIntensity, Space, Distance
    ) %>% dplyr::summarise(
        Gamma = mean(Gamma, na.rm = TRUE),
        Value = mean(Value, na.rm = TRUE)
    ), color = "red", size = 2, shape = 4
) + ggplot2::labs(
    x = "Alpha Richness",
    y = "Gamma Richness",
    subtitle = "Trajectories binned by time to equalise time steps."
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::geom_abline(slope = 1, intercept = 0)

# ggplot2::ggsave(overallalphagammabin, filename = "MNA-AlphaGammaBinned-Overview.pdf", dpi = "retina", width = 11, height = 8)

ggplot2::ggplot( 
    Diversity_AG_Binned %>% dplyr::filter(
        Environment != "Mean",
        Space != "Ring"
    ) %>% dplyr::group_by(
        Modifier, ModIntensity, Space, Distance
    ) %>% dplyr::summarise(
        Gamma = mean(Gamma, na.rm = TRUE),
        Value = mean(Value, na.rm = TRUE)
    ), 
    ggplot2::aes(
        x = Value,
        y = Gamma,
        color = interaction(Modifier, ModIntensity),
        shape = interaction(Space, Distance)
    )
) + ggplot2::geom_point(
    size = 2
) + ggplot2::labs(
    x = "Alpha Richness",
    y = "Gamma Richness",
    subtitle = "Trajectories binned by time to equalise time steps."
) + ggplot2::geom_abline(slope = 1, intercept = 0)

# ggplot2::ggsave(overallalphagammamean, filename = "MNA-AlphaGammaBinned-Means.pdf", dpi = "retina", width = 11, height = 8)

Jaccard, Space

ggplot2::ggplot(
  JaccardSpace %>% dplyr::filter(Time > 3, Space != "Ring"),
  ggplot2::aes(x = Time, y = Jaccard, 
               color = interaction(Env1, Env2)
  )
) + ggplot2::geom_line(
  alpha = 0.3
) + ggplot2::geom_line(
  data = JaccardSpace %>% dplyr::filter(Time > 3, Space != "Ring") %>% dplyr::group_by(
    Time, Distance, Modifier, ModIntensity, Space
  ) %>% dplyr::summarise(
    Jaccard = mean(Jaccard, na.rm = TRUE)
  ),
  ggplot2::aes(
    x = Time, y = Jaccard
  ),
  inherit.aes = FALSE, color = "black"
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::guides(
  color = "none"
)

Jaccard, Time

ggplot2::ggplot(
  JaccardTime %>% dplyr::filter(
    Space != "Ring"
  ) %>% dplyr::group_by(
    Environment, Distance, Space, Modifier, ModIntensity
  ) %>% dplyr::mutate(
    TimeDifference = Time2 - Time1
  ),
  ggplot2::aes(
    x = TimeDifference,
    y = Jaccard
  )
) + ggplot2::geom_bin2d(
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_fill_viridis_c(
  trans = "log"
)

ggplot2::ggplot(
  JaccardTime %>% dplyr::filter(
    Space != "Ring",
    Space == "None",
    Modifier == "Result"
  ) %>% dplyr::group_by(
    Environment, Distance, Space, Modifier, ModIntensity
  ) %>% dplyr::mutate(
    TimeDifference = Time2 - Time1
  ),
  ggplot2::aes(
    x = TimeDifference,
    y = Jaccard
  )
) + ggplot2::geom_bin2d(
) + ggplot2::facet_wrap(
    ~ Environment
) + ggplot2::scale_fill_viridis_c(
  trans = "log"
)

Species Presence

# ggplot2::ggplot(
#   SpeciesPresence %>% dplyr::filter(
#     Space != "Ring",
#     TimeFloor > 0.2 # Remove "burn-in" which has "impossibly" high presence. 
#   ),
#   ggplot2::aes(
#     x = TimeFloor,
#     y = Species
#   )
# ) + ggplot2::geom_bin2d(
#   binwidth = c(0.1, 1)
# ) + ggplot2::scale_fill_viridis_c(
# )+ ggplot2::facet_grid(
#   Modifier + ModIntensity ~ Space + Distance
# )
ggplot2::ggplot(
  SpeciesPresence %>% dplyr::filter(
    Space != "Ring"
  ) %>% dplyr::group_by(
    Modifier, ModIntensity, Space, Distance,
    Species, Time
  ) %>% dplyr::summarise(
    Count = n()
  ),
  ggplot2::aes(x = Time, y = Species, color = Count)
) + ggplot2::geom_point(
  shape = '.'
) + ggplot2::scale_color_viridis_c(
) + ggplot2::facet_grid(
  Modifier + ModIntensity ~ Space + Distance
) + ggplot2::geom_hline(
  yintercept = 34.5, color = "red"
)

ggplot2::ggplot(
  SpeciesPresence %>% dplyr::filter(
    Space != "Ring"
  ) %>% dplyr::group_by(
    Modifier, ModIntensity, Space, Distance,
    Species, Time
  ) %>% dplyr::summarise(
    Count = n()
  ),
  ggplot2::aes(x = Time, y = Species, color = Count)
) + ggplot2::geom_point(
  shape = '.'
) + ggplot2::scale_color_viridis_c(
  direction = -1
) + ggplot2::facet_grid(
  Modifier + ModIntensity ~ Space + Distance
) + ggplot2::geom_hline(
  yintercept = 34.5, color = "red"
)

# Brute forcing here because I would rather position myself to move on if possible.
# This works because all of the pools are the same (except that the original does not use factors because stringsAsFactors differed between the machines generating the systems).
# > all.equal(PoolsMats$`MNA-Ext10-PoolMats-Env10.RData`, PoolsMats$`MNA-Ext0.1-PoolMats-Env10.RData`)
# [1] TRUE
# > all.equal(PoolsMats$`MNA-Ext10-PoolMats-Env10.RData`, PoolsMats$`MNA-Arr0.1-PoolMats-Env10.RData`)
# [1] TRUE
# > all.equal(PoolsMats$`MNA-Ext10-PoolMats-Env10.RData`, PoolsMats$`MNA-Arr10-PoolMats-Env10.RData`)
# [1] TRUE
SpeciesPresence$Sizes <- PoolsMats$`MNA-FirstAttempt-PoolMats-Env10.RData`$Pool$Size[SpeciesPresence$Species]

ggplot2::ggplot(
  SpeciesPresence %>% dplyr::filter(
    Space != "Ring"
  ) %>% dplyr::group_by(
    Modifier, ModIntensity, Space, Distance,
    Species, Time, Sizes
  ) %>% dplyr::summarise(
    Count = n()
  ) %>% dplyr::arrange(
    Sizes
  ) %>% dplyr::ungroup(
  ) %>% dplyr::mutate(
    Species = factor(Species, levels = Species)
  ),
  ggplot2::aes(x = Time, y = Species, color = Count)
) + ggplot2::geom_point(
  shape = '.'
) + ggplot2::scale_color_viridis_c(
) + ggplot2::facet_grid(
  Modifier + ModIntensity ~ Space + Distance
) + ggplot2::geom_hline(
  yintercept = 34.5, color = "red"
)

Error in `levels<-`(`*tmp*`, value = as.character(levels)) : 
  factor level [2] is duplicated

Species Abundance

# ggplot2::ggplot(
#   SpeciesPresence %>% dplyr::filter(
#     Space != "Ring",
#     TimeFloor > 0.2 # Remove "burn-in" which has "impossibly" high presence. 
#   ),
#   ggplot2::aes(
#     x = TimeFloor,
#     y = Abundance
#   )
# ) + ggplot2::geom_bin2d(
# ) + ggplot2::scale_fill_viridis_c(
# ) + ggplot2::facet_grid(
#   Modifier + ModIntensity ~ Space + Distance
# ) + ggplot2::scale_y_log10(
# ) 
ggplot2::ggplot(
  SpeciesPresence %>% dplyr::filter(
    Space != "Ring"
  ) %>% dplyr::mutate(
    Abundance = round(Abundance)
  ) %>% dplyr::group_by(
    Modifier, ModIntensity, Space, Distance,
    Abundance, Time
  ) %>% dplyr::summarise(
    Count = n()
  ),
  ggplot2::aes(x = Time, y = Abundance, color = Count)
) + ggplot2::geom_point(
  shape = '.'
) + ggplot2::scale_color_viridis_c(
  trans = "log"
) + ggplot2::facet_grid(
  Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_y_log10(
)

Species Sizes

(Recall that there is a threshold above which we have consumers and below which we have producers. This will be denoted with a bright line at the threshold.)

# ggplot2::ggplot(
#   SpeciesPresence %>% dplyr::filter(
#     Space != "Ring",
#     TimeFloor > 0.2 # Remove "burn-in" which has "impossibly" high presence. 
#   ),
#   ggplot2::aes(
#     x = TimeFloor,
#     y = Sizes
#   )
# ) + ggplot2::geom_bin2d(
#   binwidth = c(0.1, 0.03)
# ) + ggplot2::scale_fill_viridis_c(
# ) + ggplot2::facet_grid(
#   Modifier + ModIntensity ~ Space + Distance
# ) + ggplot2::scale_y_log10(
# ) + ggplot2::geom_hline(
#   yintercept = 0.1, color = "red"
# )

ggplot2::ggplot(
  SpeciesPresence %>% dplyr::filter(
    Space != "Ring"
  ) %>% dplyr::group_by(
    Modifier, ModIntensity, Space, Distance,
    Sizes, Time
  ) %>% dplyr::summarise(
    Count = n()
  ),
  ggplot2::aes(x = Time, y = Sizes, color = Count)
) + ggplot2::geom_point(
  shape = '.'
) + ggplot2::scale_color_viridis_c(
) + ggplot2::facet_grid(
  Modifier + ModIntensity ~ Space + Distance
) + ggplot2::geom_hline(
  yintercept = 0.1, color = "red"
) + ggplot2::scale_y_log10(
)

Compare with the distribution of traits amongst species.

ggplot2::ggplot(
  PoolsMats$`MNA-FirstAttempt-PoolMats-Env10.RData`$Pool,
  ggplot2::aes(
    x = Size,
    fill = Type
  )
) + ggplot2::geom_histogram(#ggplot2::geom_density(
  #alpha = 0.5, adjust = 1/10
  bins = 100
) + ggplot2::scale_x_log10() + ggplot2::coord_flip()

Events

Preparation

Calculate_Events <- function(result, abundanceTime = divide_time_by, thinning = by_for_thinning) {
  abund <- result$Abundance
  abund[, -1] <- abund[, -1] > 0 # Presence-Absence
  abundDiff <- apply(abund[, -1], 2, diff) 
  # Arrival (1) Elimination (-1)
  # if arrival at time 93,
  # 92nd entry is 1.
  abundDiff <- cbind(
    abund[-1, 1] * abundanceTime, abundDiff
  )
  
  Events <- lapply(
    1:ncol(abundDiff[, -1]),
    function(i, ab, tm, sp) {
      arrivals <- which(ab[, i] == 1)
      extincts <- which(ab[, i] == -1)
      
      if (length(arrivals) + length(extincts) == 0) {
        return(NULL)
      } 
      
      data.frame(
        Times = c(tm[arrivals], tm[extincts]),
        Species = ((i - 1) %% sp) + 1, # 1:1000 -> 1:100
        Environment = ((i - 1) %/% sp) + 1,
        Type = c(rep("Arrival", length(arrivals)),
                  rep("Extinct", length(extincts))),
        stringsAsFactors = FALSE
      )
    },
    ab = abundDiff[, -1],
    tm = abundDiff[, 1],
    sp = ncol(abundDiff[, -1]) / result$NumEnvironments#,
    #ne = result$NumEnvironments
  ) %>% dplyr::bind_rows(
  ) %>% dplyr::arrange(
    Times, Species, Environment, Type
  )
  
  # Now we check to see which events are in the record.
  # Note that, due to thinning, we do have a theoretical
  # problem: an event time might be misrecorded when
  # extracted from the abundance record.
  # Hence we cannot just use filtering join operations.
  # Since we know the maximum time step size and the
  # thinning, we know that we should detect a change
  # within (thinning) * (max. time step size) units.
  maximumGap <- thinning * result$Parameters$MaximumTimeStep
  
  result$Events$Type <- as.character(result$Events$Type)
  
  # Connect all same event, even with different times.
  # Remove those that cannot be the same event.
  # Treat this as the list of Neutral Events that
  # actually happened.
  EventsOfficial <- result$Events %>% dplyr::left_join(
    Events, 
    by = c("Species", "Environment", "Type")
  ) %>% dplyr::mutate(
    # dplyr::filter( # Filtering does not work.
    # We are seeing losses of about 50% with filter.
    # "Rows in x with no match in y will have NA values in the new columns." 
    # -> Times.y will be NA.
    # Times.y is the detected time.
    # Times.x is the recorded action's time.
    # Note that we can get false readings from subtracting two almost the same 
    # numbers, so we need to appeal to machine precision. 
    # See all.equal's tolerance argument.
    `Times.y` = dplyr::case_when(
      is.na(`Times.y`) ~ as.double(NA),
      (
      `Times.y` - `Times.x` < maximumGap + sqrt(.Machine$double.eps) & 
        `Times.y` - `Times.x` >= -sqrt(.Machine$double.eps)) ~ `Times.y`,
      TRUE ~ as.double(NA)
    )
  ) %>% dplyr::group_by(
    `Times.x`, Species, Environment, Type, Success
    # Don't discard Success, others are true groups.
    # We want to preserve the first `Times.y`
    # (in case an event happens multiple times)
    # but if there are no numerics,
    # we instead want to keep one of the NAs.
  ) %>% dplyr::summarise(
    `Times.y` = if (length(na.omit(`Times.y`)) == 0) NA else min(`Times.y`, na.rm = TRUE)
  ) %>% dplyr::ungroup(
  ) %>% dplyr::rename(
    TimeImplemented = `Times.x`,
    TimeDetected = `Times.y`
  )
  
  # Events that were not detected but were successful.
  # If this is the case, something happened and was 
  # undone in the same timespan.
  # This might happen due to arrivals from adjacent 
  # patches.
  # I.e. elimination coinciding with arrival or 
  # arrival being too dissipated by dispersal.
  # This probably shouldn't happen in the disconnected system.
  EventsNotDetected <- EventsOfficial %>% dplyr::filter(
    is.na(TimeDetected), Success == TRUE
  # ) %>% dplyr::select(
  #   -`Times.y`
  # ) %>% dplyr::rename(
  #   Times = `Times.x`
  ) %>% dplyr::mutate(
    Neutral = TRUE,
    Detected = FALSE
  )
  
  # Events that were detected and were successful. 'True Positives'
  # We recorded them as having happened in Events and Abundance.
  # These are also Neutral with high probability.
  EventsDetected <- EventsOfficial %>% dplyr::filter(
    !is.na(TimeDetected), Success == TRUE
    # ) %>% dplyr::select(
    # # We keep Times.y for comparison with Events.
    #   -`Times.x`
    # ) %>% dplyr::rename(
    #   Times = `Times.y`
    ) %>% dplyr::mutate(
      Neutral = TRUE,
      Detected = TRUE
    )
  
  # Events that were not detected and were not successful. 'True Negatives'
  # Also Events that were detected and were not successful. 'False Positives'?
  # The detection must be of a different event if the event
  # we recorded was unsuccessful after all.
  EventsFailed <- EventsOfficial %>% dplyr::filter(
    #is.na(`Times.y`), 
    Success != TRUE
    # ) %>% dplyr::select(
    #   -`Times.y`
    # ) %>% dplyr::rename(
    #   Times = `Times.x`
    ) %>% dplyr::mutate(
      Neutral = TRUE,
      Detected = FALSE
    )
  
  # So we have events that were detected but not successful.
  # Such events should probably be listed twice: 
  #   once as neutral (the failure, above)
  #   once as non-neutral (the detected event, below).
  # (Note Events are from abundance and thus detected.)
  EventsNotOfficial <- Events %>% dplyr::rename(
    TimeDetected = Times
  ) %>% dplyr::anti_join(
    EventsDetected, by = c("TimeDetected", "Species", "Environment", "Type")
  )  %>% dplyr::mutate(
    Success = TRUE,
    Neutral = FALSE,
    Detected = TRUE
  ) 
  
  # The remainder of event space is events that are
  # not neutral and not successful or 
  # events that were not implemented.
  
  # Note then that, if everything went well
  stopifnot(nrow(EventsNotOfficial) + nrow(EventsDetected) == nrow(Events),
            nrow(EventsDetected) + nrow(EventsFailed) + nrow(EventsNotDetected) == nrow(result$Events))

  return(dplyr::bind_rows(
    EventsDetected,
    EventsNotDetected,
    EventsFailed,
    EventsNotOfficial
  ) %>% dplyr::mutate(
    Times = dplyr::case_when(
      !is.na(TimeImplemented) ~ TimeImplemented,
      !is.na(TimeDetected) ~ TimeDetected
    )
  ) %>% dplyr::arrange(Times, Environment, Species, Type))
}

Events <-  sapply(
  USE.NAMES = TRUE, simplify = FALSE,
  Results, function(result) {
    if (length(result) == 1 && is.na(result)) {
      # Problem case.
      return(NA)
    }
    # print(paste("Calculating", Sys.time()))
    
    # Calculate the diversity.
    # We will need to extract the system properties from
    # the file names which we carry through using sapply.
    return(Calculate_Events(result))
  }
)

Events <- lapply(1:length(Events),
                    function(i, df, nm) {
                      df[[i]] %>% mutate(
                        Simulation = nm[i]
                      )
                    }, 
                    df = Events,
                    nm = names(Events))

EventsSuccesses <- lapply(Events, function(Event) {
  Event %>% dplyr::filter(
    Success == TRUE
  ) %>% dplyr::group_by(
    Environment
  ) %>% dplyr::arrange(
    Times
  ) %>% dplyr::mutate(
    InterarrivalTime = Times - lag(Times), 
    Dynamic = !Neutral, 
    Sequence = cumsum(Neutral)
  ) %>% dplyr::group_by(
    Environment, Sequence
  ) %>% dplyr::mutate(
    EventInSequence = cumsum(Dynamic)
  )
}) 

Events <- lapply(Events, function(Event) {
  Event %>% dplyr::group_by(
    Environment
  ) %>% dplyr::arrange(
    Times
  ) %>% dplyr::mutate(
    InterarrivalTime = Times - lag(Times), 
    Dynamic = !Neutral, 
    Sequence = cumsum(Neutral)
  ) %>% dplyr::group_by(
    Environment, Sequence
  ) %>% dplyr::mutate(
    EventInSequence = cumsum(Dynamic)
  )
}) 

Events <- dplyr::left_join(
  dplyr::bind_rows(Events),
  Properties,
  by = c("Simulation" = "FullName")
)

EventsSuccesses <- dplyr::left_join(
  dplyr::bind_rows(EventsSuccesses),
  Properties,
  by = c("Simulation" = "FullName")
)

Events <- Events %>% dplyr::mutate(
  Distance = dplyr::case_when(
    is.na(Distance) ~ "1e+00",
    TRUE ~ Distance
  )
)

EventsSuccesses <- EventsSuccesses %>% dplyr::mutate(
  Distance = dplyr::case_when(
    is.na(Distance) ~ "1e+00",
    TRUE ~ Distance
  )
)

# Retrieve the Characteristic Rate used.
# Since all of the Pools and Matrices are the same, they all have the same characteristic rate.
# There is a fudge here. Looking back at the code, the first matrix was used, rather than looking over all of the matrices, but we note that the results are essentially the same.
# Perhaps something to be more careful of for the next round of simulations...
CharacteristicRate <- max(abs(eigen(PoolsMats[[1]]$InteractionMatrices$Mats[[1]])$values))

set.seed(1)
exampleData <- Events[[2]] %>% dplyr::filter(Neutral) %>% dplyr::arrange(Times) %>% pull(Times) %>% diff
fitdistrplus::descdist(exampleData, boot = 100)

censData <- ifelse(exampleData == 0, .Machine$double.eps, exampleData)
exampleFitsNotHeavy <- list(
  exp = fitdistrplus::fitdist(censData, "exp", method = "mle"),
  gamma = fitdistrplus::fitdist(censData, 
                              "gamma", method = "mle"),
  weibull = fitdistrplus::fitdist(censData, 
                              "weibull", method = "mle")
)

exampleFitsHeavy <- lapply(
    c(poweRlaw::conexp, poweRlaw::conlnorm, poweRlaw::conpl, poweRlaw::conweibull), function(f) {f$new(censData)}
)

exampleFitsHeavy <- lapply(
  exampleFitsHeavy, 
  function(d) {
    d$setXmin(poweRlaw::estimate_xmin(d, xmax = Inf))
    d
    } 
)
  
names(exampleFitsHeavy) <- c("exp", "lnorm", "pl", "weibull")

fitdistrplus::gofstat(exampleFitsNotHeavy)

quiet(exampleData_gamlss <- gamlss::fitDist(exampleData))
exampleData_gamlss

quiet(censData_gamlss <- gamlss::fitDist(censData))
exampleData_gamlss

exampleFitsHeavy

Plots

Full x Range

All Events, By Number in Sequence, Time by Environment

ggplot2::ggplot(
    Events %>% dplyr::filter(
        Space != "Ring"
    ),
    ggplot2::aes(
        x = InterarrivalTime,
        fill = factor(EventInSequence)
    )
) + ggplot2::geom_density(
  alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
  "Event #"
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10() #-> environmenteventssequence
# ggplot2::ggsave(environmenteventssequence, filename = "MNA-EventEnvironmentArrivalTimes-Sequence.pdf", dpi = "retina", width = 11, height = 8)

All Events, By Neutral/Dynamic, Time by Environment

ggplot2::ggplot(
    Events %>% dplyr::filter(
        Space != "Ring"
    ),
    ggplot2::aes(
        x = InterarrivalTime,
        fill = Neutral
    )
) + ggplot2::geom_density(
  alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10() #-> environmenteventsneutral
# ggplot2::ggsave(environmenteventsneutral, filename = "MNA-EventEnvironmentArrivalTimes-Overview.pdf", dpi = "retina", width = 11, height = 8)

All Events, By Number in Sequence, Overall Time

Events %>% dplyr::group_by(
    Simulation, Iter, Distance, 
    Modifier, ModIntensity, Environments, Space
) %>% dplyr::arrange(
    Times
) %>% dplyr::mutate(
    InterarrivalTime = Times - lag(Times), 
    Dynamic = !Neutral, 
    Sequence = cumsum(Neutral)
) %>% dplyr::group_by(
    Environment, Sequence
) %>% dplyr::mutate(
    EventInSequence = cumsum(Dynamic)
) %>% dplyr::filter(
    Space != "Ring"
) %>% ggplot2::ggplot(ggplot2::aes(
    x = InterarrivalTime,
    fill = factor(EventInSequence)
)
) + ggplot2::geom_density(
    alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
  "Event #"
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10() + ggplot2::coord_cartesian(ylim = c(0, 2)) # -> overalleventssequence
# ggplot2::ggsave(overalleventssequence, filename = "MNA-EventOverallArrivalTimes-Sequence.pdf", dpi = "retina", width = 11, height = 8)

All Events, By Neutral/Dynamics, Overall Time

Events %>% dplyr::group_by(
    Simulation, Iter, Distance, 
    Modifier, ModIntensity, Environments, Space
) %>% dplyr::arrange(
    Times
) %>% dplyr::mutate(
    InterarrivalTime = Times - lag(Times), 
    Dynamic = !Neutral, 
    Sequence = cumsum(Neutral)
) %>% dplyr::group_by(
    Environment, Sequence
) %>% dplyr::mutate(
    EventInSequence = cumsum(Dynamic)
) %>% dplyr::filter(
    Space != "Ring"
) %>% ggplot2::ggplot(ggplot2::aes(
    x = InterarrivalTime,
    fill = Neutral
)
) + ggplot2::geom_density(
    alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10() # -> overalleventsneutral
# ggplot2::ggsave(overalleventsneutral, filename = "MNA-EventOverallArrivalTimes-Overview.pdf", dpi = "retina", width = 11, height = 8)

Successful Events, By Number in Sequence, Time by Environment

ggplot2::ggplot(
    EventsSuccesses %>% dplyr::filter(
        Space != "Ring"
    ),
    ggplot2::aes(
        x = InterarrivalTime,
        fill = factor(EventInSequence)
    )
) + ggplot2::geom_density(
  alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
  "Event #"
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10(
) + ggplot2::coord_cartesian(
  ylim = c(0, 2)
) #-> environmenteventssuccesssequence
# ggplot2::ggsave(environmenteventssuccesssequence, filename = "MNA-EventSuccessEnvironmentArrivalTimes-Sequence.pdf", dpi = "retina", width = 11, height = 8)

Successful Events, By Neutral/Dynamics, Time by Environment

ggplot2::ggplot(
    EventsSuccesses %>% dplyr::filter(
        Space != "Ring"
    ),
    ggplot2::aes(
        x = InterarrivalTime,
        fill = Neutral
    )
) + ggplot2::geom_density(
  alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10() #-> environmenteventssuccessneutral
# ggplot2::ggsave(environmenteventssuccessneutral, filename = "MNA-EventSuccessEnvironmentArrivalTimes-Overview.pdf", dpi = "retina", width = 11, height = 8)

Successful Events, By Number in Sequence, Overall Time

EventsSuccesses %>% dplyr::group_by(
    Simulation, Iter, Distance, 
    Modifier, ModIntensity, Environments, Space
) %>% dplyr::arrange(
    Times
) %>% dplyr::mutate(
    InterarrivalTime = Times - lag(Times), 
    Dynamic = !Neutral, 
    Sequence = cumsum(Neutral)
) %>% dplyr::group_by(
    Environment, Sequence
) %>% dplyr::mutate(
    EventInSequence = cumsum(Dynamic)
) %>% dplyr::filter(
    Space != "Ring"
) %>% ggplot2::ggplot(ggplot2::aes(
    x = InterarrivalTime,
    fill = factor(EventInSequence)
)
) + ggplot2::geom_density(
    alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
  "Event #"
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10() + ggplot2::coord_cartesian(ylim = c(0, 2)) # -> overalleventssuccesssequence
# ggplot2::ggsave(overalleventssuccesssequence, filename = "MNA-EventSuccessOverallArrivalTimes-Sequence.pdf", dpi = "retina", width = 11, height = 8)

Successful Events, By Neutral/Dynamics, Overall Time

EventsSuccesses %>% dplyr::group_by(
    Simulation, Iter, Distance, 
    Modifier, ModIntensity, Environments, Space
) %>% dplyr::arrange(
    Times
) %>% dplyr::mutate(
    InterarrivalTime = Times - lag(Times), 
    Dynamic = !Neutral, 
    Sequence = cumsum(Neutral)
) %>% dplyr::group_by(
    Environment, Sequence
) %>% dplyr::mutate(
    EventInSequence = cumsum(Dynamic)
) %>% dplyr::filter(
    Space != "Ring"
) %>% ggplot2::ggplot(ggplot2::aes(
    x = InterarrivalTime,
    fill = Neutral
)
) + ggplot2::geom_density(
    alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10() # -> overalleventssuccessneutral
# ggplot2::ggsave(overalleventssuccessneutral, filename = "MNA-EventSuccessOverallArrivalTimes-Overview.pdf", dpi = "retina", width = 11, height = 8)

Trunc. x Range

All Events, By Number in Sequence, Time by Environment

ggplot2::ggplot(
    Events %>% dplyr::filter(
        Space != "Ring"
    ),
    ggplot2::aes(
        x = InterarrivalTime,
        fill = factor(EventInSequence)
    )
) + ggplot2::geom_density(
  alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
  "Event #"
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10(
) + ggplot2::coord_cartesian(
  xlim = c(1e-1, 1e4)
) #-> environmenteventssequence
# ggplot2::ggsave(environmenteventssequence, filename = "MNA-EventEnvironmentArrivalTimes-Sequence.pdf", dpi = "retina", width = 11, height = 8)

All Events, By Neutral/Dynamic, Time by Environment

ggplot2::ggplot(
    Events %>% dplyr::filter(
        Space != "Ring"
    ),
    ggplot2::aes(
        x = InterarrivalTime,
        fill = Neutral
    )
) + ggplot2::geom_density(
  alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10(
) + ggplot2::coord_cartesian(
  xlim = c(1e-1, 1e4)
) #-> environmenteventsneutral
# ggplot2::ggsave(environmenteventsneutral, filename = "MNA-EventEnvironmentArrivalTimes-Overview.pdf", dpi = "retina", width = 11, height = 8)

All Events, By Number in Sequence, Overall Time

Events %>% dplyr::group_by(
    Simulation, Iter, Distance, 
    Modifier, ModIntensity, Environments, Space
) %>% dplyr::arrange(
    Times
) %>% dplyr::mutate(
    InterarrivalTime = Times - lag(Times), 
    Dynamic = !Neutral, 
    Sequence = cumsum(Neutral)
) %>% dplyr::group_by(
    Environment, Sequence
) %>% dplyr::mutate(
    EventInSequence = cumsum(Dynamic)
) %>% dplyr::filter(
    Space != "Ring"
) %>% ggplot2::ggplot(ggplot2::aes(
    x = InterarrivalTime,
    fill = factor(EventInSequence)
)
) + ggplot2::geom_density(
    alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
  "Event #"
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10(
) + ggplot2::coord_cartesian(
  xlim = c(1e-1, 1e4),
  ylim = c(0, 2)
) # -> overalleventssequence
# ggplot2::ggsave(overalleventssequence, filename = "MNA-EventOverallArrivalTimes-Sequence.pdf", dpi = "retina", width = 11, height = 8)

All Events, By Neutral/Dynamics, Overall Time

Events %>% dplyr::group_by(
    Simulation, Iter, Distance, 
    Modifier, ModIntensity, Environments, Space
) %>% dplyr::arrange(
    Times
) %>% dplyr::mutate(
    InterarrivalTime = Times - lag(Times), 
    Dynamic = !Neutral, 
    Sequence = cumsum(Neutral)
) %>% dplyr::group_by(
    Environment, Sequence
) %>% dplyr::mutate(
    EventInSequence = cumsum(Dynamic)
) %>% dplyr::filter(
    Space != "Ring"
) %>% ggplot2::ggplot(ggplot2::aes(
    x = InterarrivalTime,
    fill = Neutral
)
) + ggplot2::geom_density(
    alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10(
) + ggplot2::coord_cartesian(
  xlim = c(1e-1, 1e4)
) # -> overalleventsneutral
# ggplot2::ggsave(overalleventsneutral, filename = "MNA-EventOverallArrivalTimes-Overview.pdf", dpi = "retina", width = 11, height = 8)

Successful Events, By Number in Sequence, Time by Environment

ggplot2::ggplot(
    EventsSuccesses %>% dplyr::filter(
        Space != "Ring"
    ),
    ggplot2::aes(
        x = InterarrivalTime,
        fill = factor(EventInSequence)
    )
) + ggplot2::geom_density(
  alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
  "Event #"
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10(
) + ggplot2::coord_cartesian(
  xlim = c(1e-1, 1e4),
  ylim = c(0, 2)
) #-> environmenteventssuccesssequence
# ggplot2::ggsave(environmenteventssuccesssequence, filename = "MNA-EventSuccessEnvironmentArrivalTimes-Sequence.pdf", dpi = "retina", width = 11, height = 8)

Successful Events, By Neutral/Dynamics, Time by Environment

ggplot2::ggplot(
    EventsSuccesses %>% dplyr::filter(
        Space != "Ring"
    ),
    ggplot2::aes(
        x = InterarrivalTime,
        fill = Neutral
    )
) + ggplot2::geom_density(
  alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10(
) + ggplot2::coord_cartesian(
  xlim = c(1e-1, 1e4)
) #-> environmenteventssuccessneutral
# ggplot2::ggsave(environmenteventssuccessneutral, filename = "MNA-EventSuccessEnvironmentArrivalTimes-Overview.pdf", dpi = "retina", width = 11, height = 8)

Successful Events, By Number in Sequence, Overall Time

EventsSuccesses %>% dplyr::group_by(
    Simulation, Iter, Distance, 
    Modifier, ModIntensity, Environments, Space
) %>% dplyr::arrange(
    Times
) %>% dplyr::mutate(
    InterarrivalTime = Times - lag(Times), 
    Dynamic = !Neutral, 
    Sequence = cumsum(Neutral)
) %>% dplyr::group_by(
    Environment, Sequence
) %>% dplyr::mutate(
    EventInSequence = cumsum(Dynamic)
) %>% dplyr::filter(
    Space != "Ring"
) %>% ggplot2::ggplot(ggplot2::aes(
    x = InterarrivalTime,
    fill = factor(EventInSequence)
)
) + ggplot2::geom_density(
    alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
  "Event #"
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10(
) + ggplot2::coord_cartesian(
  xlim = c(1e-1, 1e4),
  ylim = c(0, 2)
) # -> overalleventssuccesssequence
# ggplot2::ggsave(overalleventssuccesssequence, filename = "MNA-EventSuccessOverallArrivalTimes-Sequence.pdf", dpi = "retina", width = 11, height = 8)

Successful Events, By Neutral/Dynamics, Overall Time

EventsSuccesses %>% dplyr::group_by(
    Simulation, Iter, Distance, 
    Modifier, ModIntensity, Environments, Space
) %>% dplyr::arrange(
    Times
) %>% dplyr::mutate(
    InterarrivalTime = Times - lag(Times), 
    Dynamic = !Neutral, 
    Sequence = cumsum(Neutral)
) %>% dplyr::group_by(
    Environment, Sequence
) %>% dplyr::mutate(
    EventInSequence = cumsum(Dynamic)
) %>% dplyr::filter(
    Space != "Ring"
) %>% ggplot2::ggplot(ggplot2::aes(
    x = InterarrivalTime,
    fill = Neutral
)
) + ggplot2::geom_density(
    alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
) + ggplot2::facet_grid(
    Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10(
) + ggplot2::coord_cartesian(
  xlim = c(1e-1, 1e4)
) # -> overalleventssuccessneutral
# ggplot2::ggsave(overalleventssuccessneutral, filename = "MNA-EventSuccessOverallArrivalTimes-Overview.pdf", dpi = "retina", width = 11, height = 8)